﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary
{
    /// <summary>
    /// This class represents a m ultidimensional Taylor series.
    /// </summary>
    [Serializable]
    public class MultidimensionalTaylorSeries : IMultidimensionalTaylorSeries
    {
        /// <summary>
        /// Internal value for the function evaluation at position f(a).
        /// </summary>
        private readonly double fa;

        /// <summary>
        /// The point around the series approximates the function.
        /// </summary>
        private GeneralVector a;

        /// <summary>
        /// This field holds the jacobian vector.
        /// </summary>
        private GeneralVector jacobianVector;

        /// <summary>
        /// This field holds the hessian vector.
        /// </summary>
        private Matrix hessianMatrix;

        /// <summary>
        /// The function, which the series approximates.
        /// </summary>
        private HardMultivariateRealFunction function;

        /// <summary>
        /// Initializes a new instance of the <see cref="MultidimensionalTaylorSeries"/> class.
        /// </summary>
        /// <param name="a">The point around the series approximates the function.</param>
        /// <param name="jacobianVector">The jacobian vector of the taylor series.</param>
        /// <param name="hessianMatrix">The hessian matrix of the taylor series.</param>
        /// <param name="function">The function, which the series approximates.</param>
        public MultidimensionalTaylorSeries(GeneralVector a, GeneralVector jacobianVector, Matrix hessianMatrix,
                                            HardMultivariateRealFunction function)
        {
            if (function == (HardMultivariateRealFunction) null)
            {
                throw new ArgumentNullException("function");
            }

            this.a = a;
            this.function = function;
            this.fa = function.SolveAt(a);
            this.hessianMatrix = hessianMatrix;
            this.jacobianVector = jacobianVector;
        }

        /// <summary>
        /// Gets or sets the point around the series approximates the function.
        /// </summary>
        /// <value>The point around the series approximates the function.</value>
        public GeneralVector A
        {
            get { return this.a; }
            set { this.a = value; }
        }

        /// <summary>
        /// Gets or sets the jacobian vector of the taylor series.
        /// </summary>
        /// <value>The jacobian vector of the taylor series.</value>
        public GeneralVector JacobianVector
        {
            get { return this.jacobianVector; }
            set { this.jacobianVector = value; }
        }

        /// <summary>
        /// Gets or sets the hessian matrix of the taylor series.
        /// </summary>
        /// <value>The hessian matrix of the taylor series.</value>
        public Matrix HessianMatrix
        {
            get { return this.hessianMatrix; }
            set { this.hessianMatrix = value; }
        }

        /// <summary>
        /// Gets or sets the function, which the series approximates.
        /// </summary>
        /// <value>The function, which the series approximates.</value>
        public HardMultivariateRealFunction Function
        {
            get { return this.function; }
            set { this.function = value; }
        }

        /// <summary>
        /// Solves the specified Taylor series at the specified position x.
        /// </summary>
        /// <param name="x">The specified position x.</param>
        /// <returns>The series value at the specified position x.</returns>
        public double SolveAt(GeneralVector x)
        {
            if (this.a.Count != this.jacobianVector.Count)
            {
                throw new ArgumentException
                    ("The number of components of vector a and the first partial derivative have to be even");
            }

            GeneralVector xa = x - a;

            return this.fa + xa*this.jacobianVector + 0.5*xa*(this.hessianMatrix*xa);
        }
    }
}